Hadronization Approach for a Quark-Gluon Plasma Formed in Relativistic Heavy Ion 

Collisions 
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A transport model is developed to describe hadron emission from a strongly coupled quark-gluon 
plasma formed in relativistic heavy ion collisions. The quark-gluon plasma is controlled by ideal 
hydrodynamics, and the hadron motion is characterized by a transport equation with loss and gain 
terms. The two sets of equations are coupled to each other, and the hadronization hypersurface is 
determined by both the hydrodynamic evolution and the hadron emission. The model is applied to 
calculate the transverse momentum distributions of mesons and baryons, and most of the results 
agree well with the experimental data at RHIC. 
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I. INTRODUCTION 

The hot and dense matter created in relativistic heavy 
ion collisions provides the condition to investigate the 
new state of matter of QCD, the so-called quark-gluon 
plasma(QGP) 0,0,13,11,8. 

The space-time evolution of 
QGP and its hadronization process are still open ques- 
tions in the research of this area. Hydrodynamics [y, [D, [8| 
based on the assumption of local thermal and chemical 
equilibrium, models of quark coalescence or recombina- 
tion [IE 111 El 13 El, and fragmentation [H,E1 are 
often used to describe the evolution and hadronization of 
QGP formed in heavy ion collisions at RHIC. However, 
we do not know a priori if the colliding system can reach 
equilibrium or not. In fact, because of the very small size 
and very short lifetime of the collision zone [IT], EM Ell , 
the highly excited system may spend a considerable frac- 
tion of the lifetime in a non-equilibrium state (20| . On 
the other hand, even if the system reaches equilibrium in 
the early quark-gluon stage, the hadrons formed in the 
later stage will be in non-equilibrium state, since the par- 
ticle density in hadron phase is much smaller than that 
in the parton stage. 

Evaporation models [2l|, [H developed many years 
ago are widely used to describe the hadronization of 
QGP. In some of the models, the emission rate is cal- 
culated from a simple assumption of balance between 
the particle emission and absorption in QGP and in 
hadron gas at a critical temperature T c . Under this as- 
sumption, the particle emission distribution is just the 
corresponding ideal thermal distribution of hadron gas, 
/(P) = f T iP)- 

However, considering the feedback effect of hadron 
evaporation on QGP, the space-time evolution of QGP 
and hadron evaporation are coupled to each other and the 
system is governed by both the hydrodynamic expansion 
and the evaporation. On the other hand, when the ex- 
pansion of the system is fast, the fluid will be more likely 
to split into relatively small bubbles during the hadron 
evaporation. In this case, the system is a mixture of QGP 
bubbles and physical vacuum, and the energy density of 



the system is lower than that of a unitary QGP. These 
bubbles may provide a mechanism to reduce the viscosity 
of the system and make it more accessible for the ideal 
hydrodynamics to work at RHIC 24]. Thirdly, from the 
analysis of collective flow measured at RHIC (25|, the 
QGP formed at RHIC is a strongly coupled quark-gluon 
plasma (sQGP) [26|, [27], [H, |2^| where there exist not only 
quarks and gluons but also some of their resonant states 
like qq and qg. 

In this paper, we consider an effective hadronization 
model for the production of soft particles in relativistic 
heavy ion collisions. In the model, the QGP evolution 
is described by ideal hydrodynamics together with the 
equation of state including quarks, gluons and their res- 
onant states. At the critical temperature, the expansion 
does not change the temperature of the system, and the 
decrease of the energy density is due to the splitting of the 
QGP fluid. When the energy density reaches the value for 
hadronization, which is a free parameter in our calcula- 
tion, the evaporation starts and the hadron distributions 
are controlled by transport equations. The evaporation 
hypersurface is characterized by both the hydrodynamics 
and the evaporation. 

The paper is organized as follows. In Section [TT] we 
discuss hadron production and transport and derive an 
approximate emission distribution. In Section lHTl we con- 
sider the hydrodynamics for the evolution of QGP with 
possible splitting into droplets and the effective evapora- 
tion boundary condition. We apply our model to hadron 
spectra in heavy ion collisions at RHIC energy in Section 
IIVI We summarize in Section [V] 



II. HADRON EMISSION 

We investigate in this section the transport of hadrons 
emitted from QGP. The classical transport equation for 
the evolution of hadron phase space distribution f(t, r, p) 
can be written as 



df 
dt 



vVf = Q + a-0f, 



(1) 



2 



where the second term on the left hand side is the free 
streaming part, and on the right hand side the lose and 
gain terms (3(p)f and a(p) indicate hadron absorption 
and production in QGP, and the elastic term Q(p; /) [30j 
controls the thermalization of the hadrons. 

Under the assumption of quasi-static condition, the 
time evolution f3j] during the transport can be tem- 
porarily removed. If we only focus on a one dimensional 
hadron transport in the region x > and treat the elastic 
scattering in relaxation time approximation, the hadron 
transport can be greatly simplified, 



v cos 9 



df 
dx 



f-f 1 



(2) 



where v = p/E is the hadron velocity, r the relaxation 
time parameter, and f T (p) the thermal hadron distribu- 
tion. In principle, r should depend on the hadron mo- 
mentum and the type of hadrons. For simplicity, we take 
it as a universal constant in our global model. Consid- 
ering the finite volumes of hadrons and resonances, the 
phase space distribution will be suppressed. An effective 
way to include the volume effect in the thermal distri- 
bution is to multiply f T by a factor 1 — y 0) where y is 
the ratio of the occupation volume of all the excluded 
hadrons to the whole volume of the system [3^, [33[ . In 
the following calculation we will always consider this vol- 
ume effect for f T . 

The transport equation ^ can be solved with the an- 
alytic result 

f{x,e,p) = fM + [f%6 i p)-f 00 (p)}e-^r* (3) 



with fcaip) defined as 



fooip) 



f T (p) + a(p)T 
p(p)r + 1 



(4) 



It is easy to see that the hadrons reach chemical equilib- 
rium in the limit of vanishing inelastic interactions with 
a = f3 — and thermal equilibrium in the limit of in- 
finitely strong elastic interaction with r = 0. In the both 
cases, the distribution at x — > oo is reduced to the ther- 
mal one, 



fM = f T (p)- 



(•5) 



Taking into account the boundary condition at x = 0, 

/(O,0,p)=O, (cos6»>0), (6) 

which means no hadrons moving from the vacuum in the 
region of x < to the QGP phase in the region of x > 0, 
and the constraint 

0</(*,M</cc(p), (7) 
the hadron distribution can be expressed as 

f(x,9,p) = JUp) [l - etcos^e-^H , (8) 



where Q(x) is a step function. 

It can be seen from the above discussion that the in- 
fluence of the transport process is mainly reflected in the 
region near the QGP surface at x = 0. Most of the emit- 
ted hadrons are created in this region. Since the energy 
density in the inner part of the QGP fluid is higher than 
the value on the surface, the hadronization of the QGP 
fluid can be effectively regarded as a surface evaporation. 

In the above derivation we did not consider the QGP 
surface shift and the breaking of the fluid in the process 
of hadron emission. A simple way to take the surface 
shift and fluid splitting into account is to introduce two 
parameters v s and 77 in the model. The former is the 
shrinking velocity of the surface, and the latter is the 
ratio of the occupation volume of all the QGP droplets 
to the whole volume of the system. These two param- 
eters should be considered simultaneously in the trans- 
port equation for hadron emission and the hydrodynamic 
equation for QGP evolution. With these two parameters, 
the hadron transport equation can be effectively written 
as 



(vcosi 



Vs) d~x =11 



f-f 



+ a-0f 



(9) 



with the solution 

f(x,e,p) = /oo(p) x 

1 — <d(vcos8 — v s )e~ 



(10) 



With the hadron distribution, the energy loss of the 
QGP fluid in the rest frame of the boundary can be ex- 
pressed as 



dE 
dtdS 



E 



fooiPi) ( v * cos 6» - v s ) E^p,, (11) 



where the integration over the angle 9 is restricted by the 
condition Vi cos 6 — v s < 0, Ei is the energy of a hadron 
of type i, and the summation runs over all the hadron 
types. Similarly, the momentum loss of the QGP fluid 
can be written as 



E / ftoiPi) (viCos6 -v s )pcos6d 3 pi. (12) 



dP 
dtdS 



If hydrodynamic expansion is neglected, the equation 



dE 
dtdS 



(13) 



defines a minimal shrinking velocity of the fluid, where 
e is the average energy density at the boundary. When 
e is small enough to make v s = 1, the local association 
between different parts of the fluid will completely vanish, 
and the fluid must decouple. Therefore, the condition 
for the fluid description of the system is that the energy 
density must satisfy 



e > 



dE 



dtdS 



(14) 
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The question left for the emission distribution of 
hadrons is the determination of the loss and gain terms (3 
and a. We first calculate the absorption rate of hadrons 
by the medium. Considering that tens of kinds of hadrons 
will be discussed, the calculation of their different absorp- 
tions would be quite complicated. For simplicity, the in- 
elastic absorption rate of a hadron by partons of type j 
is described as 



Pj(p) = J Sy(q)|v- v q \a j {p,q)dq 



(15) 



A f 2 / N f s( V + Vq) 3 - \V-V q 

J 6vv q 

where v q = q/E denotes the parton velocity, g 3 is 
the parton thermal distribution at the critical temper- 
ature T c determined by the hydrodynamic equation, and 
erj(p, q) is the effective absorption cross section. Replac- 
ing <7j (p, q) by its average value (Xj, we obtain the total 
absorption rate 



with 



CM = 



vgj(q)q 2 + ^9j(q)jp 



(16) 



9j{q) 



-g 3 (q)qE 3 



dq 
dq. (17) 



Since the diquark dynamics is not yet very clear, to sim- 
plify the calculation we take the assumption that the 
baryon cross section is 1.5 times the meson cross section, 

(<T b ) = 1.5(<7 m ). 

For high momentum hadrons passing cold matter, we 
have 

f°° 1 f°° q 4 - 

Cjip) ~ v / 9j{q)q 2 dq+ — / gj(q)—^dq 
Jo Ju Jo a i 



•"3 
ill 



6V J 



(18) 



where, (in) is a constant at given temperature and chem- 
ical potential, and rij is the parton density. In this case, 
the total rate can be approximately expressed as 



(19) 



where L = 1/ '^2 i n j&j is the inelastic mean free path. 
For soft hadrons in hot medium, we have 

f°° q 3 v 2 f°° 

CjiP) ~ / 9j(q)-^rdq+— g j {q)qE j dq 
Jo 6 Jo 



j 

4tt 



3 Vi 



(20) 



Taking (vj 



1 for gluons and light quarks and 



(vj) w 0.96 and (— ) « 1.05 for strange quarks with mass 
m s — 95 MeV, the total rate becomes approximately 



Pip) 



l + v 2 /3 



(21) 



There are many mechanisms to describe hadron pro- 
duction in heavy ion collisions, such as two- or three- 
quark recombination for mesons or 

baryons[H[Il[I3. 

In our model, there are not only quarks but also diquarks 

in the QGP. We assume that, mesons and intermediate 
baryons are respectively made of two quarks and one 
quark and one diquark, and final state baryons are con- 
sidered as combinations of different intermediate baryon 
states [26] , see Appendix 1X1 

The produced hadron density in coordinate and time 
space through combination of two constituents j and k 
at a fixed channel / can be defined as 1341 



ni 



J h ■ - vk\gj{pj)gk(pk)(r{pj,pk,p)d 3 p j d 3 pk, (22) 



where gj and gk are the constituent distributions, and a 
is the cross section, 



<r(Pj,Pk,P) = (2tt) 4 



\Mrf 



4\/ (Pj ■ Pk) 



2 9 



(23) 



xS 4 ( Pj + P k-p)S(p 2 - m 2 )Q{E)-^ 

with the scattering matrix element Mj. 

For baryon production, since a diquark is much heavy 
than the corresponding quark, the hadron density can be 
approximately expressed as [26|, [35| , 



m = 2^ 3 |M / | 2 / gjiE^gkiE^dEjdEk, (24) 



where the integrated region is restricted by the con- 
straints 



rrij < Ej < co, 

{E 3 +Ek){m 2 Ek+mlE,j) 



(25) 



EjEk > 

E c = 



El 



— — \J [m 2 — (rrij + mk) 2 ] [m 2 — (mj — m^) 2 



Since we have assumed thermalization of quarks, glu- 
ons and diquarks, their distribution functions can be 
written as 



9j(Pj) = 9j(Ej) 



(26) 



where ujj = 3Jj(Jj + 1) is the degeneracy factor of spin 
and color degrees of freedom, and k is the reduction fac- 
tor resulted from the volume effect of the bound-state or 
resonance-state components [331 ] . 



K=[(l-y ) + C t y ] 2 - l (l-y o y, 
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with I being the number of bound states or resonant 
states and C t the "transparency" of bags. For simplicity, 
we set C t = 1. 

Making summation over all the possible channels, we 
obtain the final hadron density 



(27) 



and qg states, their masses in the vacuum [33[ are taken 

580 MeV in 0+ channel, 
820 MeV, m ss = 905 MeV in 

1 + channel, and m qg = 610 MeV, m sg = 715 MeV in \ 
channel. These values are a bit smaller than the estima- 
tion in constituent quark model (39[ . 



as m qq — 485 MeV, m qs 
m qq = 725 MeV, m qs 



where Cj is the weight factor of channel / and is listed 
in Appendix [Al for different channels. 

To have detailed information on the produced hadrons, 
we need the hadron momentum distribution 



dni 



7T 2 \M, 



A'Kp 2 dp 2 pE 
where F is defined as 



^- J F(E k ) 9j (E-E k )g k (E k )dE k , 



(28) 



EE k {m< 



{E 2 m 2 



- Elm?) 



mEr 



in r. 

and the integrated region is restricted by the constraint 
E max = Max{m k ,E-}, 

Ernin = Min{E - TTlj , E + } , 

E (m 2 - m 2 + mj) ± 2mpE c 

E± = 7T^ ■ \ M ) 

2m z 

With the known produced hadron density, the gain 
term a in the transport equation ([2]) can be written as 



a(p) = Y,Ci 



dnj 



(ATrp 2 )dp 



(31) 



Note that in the calculation of lose and gain terms 
above, we considered only hadronization and hadron in- 
teractions with the QGP constituents. When the hadron 
decay in the transport process is included, it will cer- 
tainly change the loss and gain terms. Its contribution 
to the production rate a will be effectively discussed in 
Section HVl and the change in the absorption rate (3 can 
be described by a shift 36] 



7 



(32) 



where T is the hadron decay width and 7 is the Lorentz 
factor of the hadron. 

In principle, all the particle masses are temperature 
and chemical potential dependent [13, l38| . For quarks 
and gluons at the hadronization, we can simply take 
m g = 0, m u = 3 MeV, to^ = 6 MeV and m s = 95 
MeV [36|, but for light hadrons, their masses are in prin- 
ciple controlled by chiral symmetry and therefore deviate 
from their vacuum values. To simplify calculations, we 
will still use their masses and decay widths in the vac- 
uum and neglect the medium effect. As for the diquark 



III. QGP EVOLUTION 

We now turn to the space-time evolution of QGP. For 
ideal QGP, it can be described by the hydrodynamic 
equations, 



dTu* = 0, d»jl 



0. 



(33) 



where T M „ = (e+P^^u^ — Pg^ is the energy-momentum 
tensor of the fluid, j* = riiU^ refers to the conserved 
charge current like baryon current and strange current, 
and e, P, rii and are respectively the energy density, 
pressure, conserved charge density and four velocity of 
the fluid. 

To close the hydrodynamic equations, we need the 
equation of state of the fluid, e(P). Since the parton fluid 
formed at RfflC energy is a strongly coupled QGP, we 
include in the equation of state the it, d, s quarks, gluons, 
and all possible qq and qg resonant states. To simplify 
the calculation, we assume the relation e = 3P in the 
QGP stage. 

From the hydrodynamic equation (|33[) and the equa- 
tion of state, we can obtain the energy density e as a func- 
tion of space and time. When e reaches e c corresponding 
to the critical temperature T c for deconfinement, the ex- 
pansion of the system will no longer change the tempera- 
ture, and the decrease of e is assumed due to the splitting 
of the QGP fluid. Therefore, for the expansion process 
at fixed temperature T c , the ratio of the energy densi- 
ties e/e c is always less than 1. Finally, when the energy 
density drops down to the value £(, for the hadronization, 
which is a free parameter in our calculation, the hadron 
emission starts and the hydrodynamic evolution of the 
system becomes 



8» (T M „ - TV) = 0, 



-3 U =0, 



(34) 



where T ^ and j' are the changes in energy-momentum 
tensor and conserved charge current due to the hadron 
emission, and they are determined by the energy loss (jll[) 
and momentum loss (|12j) . In this way, the evolution of 
the decoupling hypersurface of the fluid is determined by 
both the hydrodynamics and the evaporation, similar to 
the method in FCT calculations [40l fill]. 

We discuss now, as an example, the spherical expan- 
sion of hydrodynamics to demonstrate the influence of 
evaporation on the hadronization boundary. The 1 + 1 



5 



hydrodynamics of QGP can be written as [4 
2 



at 
dv 



Of 



i-v 2 

1 

Z-v 2 



df dv , . v 

v ir + 2 ir + 2 ( N - 1 )- 

or or r 
3(1 - v 2 ) 2 <9g dv 
4 <9r 6V 



-(N-l) 



v 2 (l - v 2 ) 



dr 7 dr 



dv 

+ n 7 — + (N — l)n. 



(35) 



for the spherical hydrodynamic velocity v, the scaled 
energy-momentum quantity £ = In [(e + P)/(e& + Pb)], 
and the conserved charge density n 7 = 7J_b7, where 7 is 
the Lorentz factor of the fluid, rig is the baryon density, 
and TV = 3 is the symmetric dimension. 

Similarly, we can derive the simplified 1+1 hydrody- 
namic equation including the feedback from the hadron 
evaporation. When e calculated from the hydrodynamics 
with evaporation is less than e&, the hadronization hyper- 
surface is determined by ([35f . Otherwise, the hydrody- 
namics with evaporation itself gives the hypersurface. 

The radius R of the hadronization spherical surface, 
the hydrodynamic expansion velocity v on the surface, 
and the bulk energy of the fluid are shown as functions 
of time in FigfT] for different values of shrinking veloc- 
ity v s which characterizes the effect of hadron emission 
on the shift of hydronization hypersurface. In the cal- 
culation, we have chosen the initial energy density dis- 
tribution as a Gaussian distribution with the central 
value e(r = 0) = 9e;, and the initial expansion veloc- 
ity u(r) = 7W = 0.05r fm^ 1 . For v s = which means 
no hadron emission, the radius and the expansion veloc- 
ity are fully determined by the hydrodynamic equations 
(|35[) . With increasing hadron emission, the shrinking ve- 
locity increases, and the hadronization radius, the expan- 
sion velocity and the bulk energy deviate from the pure 
hydrodynamic solution more and more strongly. Since 
the hydrodynamic expansion grows with increasing time, 
the quantities are always controlled by the pure hydro- 
dynamics and the emission effect can be neglected in the 
later stage of the fluid. 

To study the space-time evolution of relativistic heavy 
ion collisions at RHIC energy, it is better to consider the 
cylindrical expansion of hydrodynamics. As is usually 
done, the longitudinal expansion can be decoupled from 
the transverse expansion and described in the Bjorken 
boost invariant scenario [43^ . The transverse expansion 
is described by equation (|35f with N = 2. As an ex- 
ample, we consider a central Au-Au collision at colliding 
energy ^/s — 200 GeV. The initial condition of the col- 
lision, which is the same as in [f| [l{| , is determined by 
the nuclear geometry [lij and the corresponding nucleon- 
nucleon interaction. The formation time of the locally 
thermalized fluid is set to be r = 0.6 fm/c. 

The hadronization radius R in the transverse plane and 
the corresponding expansion velocity v at R are shown 
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FIG. 1: The hadronization radius R, the expansion velocity v 
at the radius, and the bulk energy E scaled by the hadroniza- 
tion energy density £5 of the QGP fluid as functions of time 
for several values of shrinking velocity v 3 , calculated in the 
frame of spherical expansion. 



in Figf2] The hadronization energy density is chosen as 
e b = 0.21 GeV/fm 3 by fitting the RHIC data, and the 
shrinking velocity v s is the maximum of its value solved 
from the self-consistent equations (fTTj) and (fT3"|) for the 
emission energy and its value determined from the QGP 
hydrodynamic equation with splitting. Different from 
the spherical expansion where the hadron emission plays 
an essential role in determining the hadronization radius, 
the evaporation in the frame of cylindrical expansion does 
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not change the pure hydrodynamic result considerably, 
except in the very beginning period. This behavior is 
due to the rapid expansion in the Bjorken scenario. 





FIG. 2: The hadronization radius R and the corresponding 
expansion velocity v of the QGP fluid as functions of time in a 
central Au-Au collision with colliding energy y/s — 200 GeV 
per pair of nucleons, calculated in the frame of cylindrical 
expansion. 



IV. HADRON PRODUCTION AT RHIC 

Now it is the time for us to apply our model with hy- 
drodynamics for QGP evolution and transport for hadron 
emission to relativistic heavy ion collisions. We will cal- 
culate the spectra of soft hadrons produced in central 
Au-Au collisions at colliding energy ^snn = 200 GeV 
and compare the results with experimental data at RHIC 

HE EE m, HI Hi M M, Hi M 0, HE M, i3 ■ 

Let's first list the free parameters in the model. We 
choose the critical temperature T c — 166 MeV which 
leads to the critical energy density e c = 1.91 GeV/fm 3 . 
Since resonant states of qq and qg are taken into account 
in our equation of state for the strongly coupled QGP, the 
energy density is larger than the Stefan-Boltzmann limit 
for the ideal system of quarks and gluons, e c — 1.2Aesb- 
The average cross section for meson-quark inelastic scat- 



tering is taken to be (a m ) = 7 mb. The ratio y to de- 
scribe the finite volume effect of hadrons is taken as 0.83. 
From our numerical calculation, the obtained hadron 
spectra are not sensitive to the value of y in a wide re- 
gion. To minimize the free parameters, we assume that 
the hadrons are not thermalized at all, namely we take 
the relaxation time r = oo for all hadrons. In this case, 
the hadron distribution is controlled only by the emission 
rate and medium absorption, 



/oo(p) 



0{p)' 



(36) 



Further more, by fitting the PHENIX data [45( of pro- 
tons and pions, we determine the energy density of the 
QGP fluid e& — 0.21 GeV/fm 3 at hadronization which 
corresponds to the ratio r\ = 6b/e c = 0.11, the chemical 
potential pb = 27 MeV at hadronization, and the meson- 
quark and baryon-quark coupling constants gj^ — 5.11 
and g\ = 34.24. For strange quark s, its chemical po- 
tential is fi s = fis/3 — [is where /is is the strange chem- 
ical potential. With the above parameters, the evolution 
duration of the QGP fluid, calculated from the hydrody- 
namic equation, is about 8.4 fm/c. 

The invariant yields of hadrons with type i can be cal- 
culated from the phase space integration of the hadron 
distribution 



d 2 N, 



2npTdpTdy 2-KpxdpT 



E / /i(q')(^c°s0-< 



xF 1 [ ;(( i ,p)d z xdt'd i (i, 



(37) 



where Fp stands for the decay from hadrons of type j to 
the hadrons of type i, q' is the hadron momentum in the 
local rest frame of the hadronization hypersurface and is 
related to the momentum q by a Lorentz transformation. 
It is necessary to note that, the decay contribution to 
the gain term a, mentioned in the last section, can be 
effectively described by multiplying the decay term Fp 
here by a factor 



a ; 



with 



1 + £.7+1,3 + £j + lj"G+2j + l H h 0+1 ,j ■ ■ ■ 

(38) 







+1,3 



03+1 + r j + l/ ' ilj+lti) ' 



(39) 



where we have considered the assumption of r = oo. 

The phase space integration is done by a Monte-Carlo 
simulation, which partially trace the directly emitted 
hadrons. By counting enough events, the statistical error 
for dN/dy caused by the random simulation is less than 
0.5% — 1% for different hadrons. The statistical error for 
Pt distribution is, however, relatively larger, especially 
in very low px and high px regions. The reason is the 
denominator pt in (|37p for low px hadrons and low abun- 
dance for high pt hadrons. In the region of px — 0.3 — 4 
GeV, the statistical error is expected to be smaller than 
5% for most hadrons. 
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The light hadron yield is dominated by decay product. 
For example, the number of direct pions is less than 5% 
of the total final state pion abundance. Since the con- 
vergence of heavy resonance decay contribution is slow, 
many resonances need to be considered in our calcula- 
tions. The selection of hadron types is subject to their 
relative production rates and branch ratios, compared 
with their daughter particles counted in our calculations. 
Mesons and baryons included in our model are listed in 
Table HI and [TT1 Although some of + particles, such as 
/o(600), / (980) and a (980) are often interpreted as non- 
qq states, they are still treated as scalar mesons in our 
calculations. 
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FIG. 3: The transverse momentum distributions of 7T . The 
abundances are inclusive, and no feed-down corrections are 
applied. The data are from PHENIX [U and PHOBOS [I| 
Collaborations. 



TABLE I: Mesons selected in the calculations. 
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A 


A 


E 




Q 


Ng3s 


A1232 


Ain6 


E1192 


Sl317 


^1672 


A r i440 


A1600 


A1405 


El385 


Sl530 




A r i520 


A1620 


A1520 


El660 


Hl690 




^1535 




A16OO 


£l670 


Hl820 








Al670 


£l750 










Al690 


El 775 










Al830 









TABLE II: Baryons selected in the calculations. 
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The hadron spectra and relative ratios are calculated 
with the invariant momentum distribution (I37p at y = 0. 
Since the data for exclusive protons and inclusive pions 
from PHENIX [45| are used to fit the parameters, the 
agreement between our model calculation and the data 
for these two kinds of particles, shown in Figs G2 2] and 
looks perfect. 

As only soft particles are emitted from the QGP fluid 
itself, the hadron spectra in high transverse momentum 
region are expected not to be reproduced correctly in our 
model. For these high pr particles, we should consider 
the jet motion in QGP and its fragmentation mechanism. 
The effect of jets or shower partons in high pt region is 
more significant for pions and protons, compared with 
those heavy resonances, especially multi-strange heavy 
resonances. The big difference between our model calcu- 
lation and STAR data [13] for high pt pions and protons 
is shown clearly in Figs. [6] and [71 

The transverse momentum distributions for charged 
kaons are shown in Fig. [5] Our model result looks higher 
than the experimental data. The reason is probably the 



FIG. 4: The transverse momentum distributions of protons 
and anti-protons. The feed-down corrections from A are ap- 
plied. The data are from PHENIX Collaboration 



lack of heavy resonances in our hadron list or the light 
strange quark with m s — 95 MeV 3G] . A heavier strange 
quark can reduce the production of strange hadrons like 
kaons, E*(1385) and E*(1530). For m s = 170 MeV, the 
suppression for kaons is about 5%. Another possible rea- 
son for the overestimation of kaons is the strange hadron 
thcrmalization. If we take the relaxation time r ~ 1 fm, 
the model calculation agrees well with the data. 

The A and S~ spectra are shown in Figs. [9] and 
mil The model describes well the S data, but the re- 
sult deviates from the data for A and A remarkably in 
low pt region, especially for the positive A. Compared 
with the experimental data, a general behavior of our 
model calculation is the enhancement of strange mesons 
and suppression of strange baryons. This is probably 
due to the exclusion of the hadron rescattering, such as 
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FIG. 5: The p/n + and p/ir ratios as functions of transverse 
momentum. The data are from PHENIX Collaboration 14511. 
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FIG. 6: The transverse momentum distributions of proton 
and 7r + . The distribution for proton is inclusive, and the 
decay contribution from A and Kg to p ions is removed. The 
data are from STAR Collaboration [47| . 



FIG. 7: The p/tt + ratio as a function of transverse momen- 
tum. The distribution for proton is inclusive, and the decay 
contribution from A and K s to pions is removed. The data 
are from STAR Collaboration [47[. 
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FIG. 8: The transverse momentum distributions of K . The 
abundances are inclusive and no feed-down corrections are 
applied. The data are from PHENIX Collaboration [4jJ. 



N + K *=> E(A) + tt or N + K —> E*(A*) + AT, or the 
treatment of diquarks. In our current frame, diquarks 
are thermalized like quarks and gluons and satisfy the 
Bose-Einstein distribution. Since diquarks are so heavy, 
it is more reasonable to describe their phase-space distri- 
bution by a transport equation, like the one for hadrons. 

The transverse momentum distributions for A, K + and 
Kg and the ratio A/ Kg are shown in Figs. [TT1 and [T2"l 
Again we see clearly the theoretical enhancement for K 
and suppression for A. From the comparison with the 
STAR data, the inclusive A distribution without the feed- 
down correction seems much better than the one shown 
in FigEl Different from the theoretical soft p/ir ratio 
shown in Fig[7| where the ratio never decreases, the soft 
A/K ratio begins to drop down at about pt ~ 4 GeV. 

The if* (892) yield is shown in Fig. [T3] as a function 



of m,T — ttik*i where m-r — \JVt ^ m k* ^ s ^ ne P ar ticle 
transverse mass. From the comparison with the STAR 
data, the decay correction from the loss term (|3"2")) plays 
an important role. As the decay products may be ab- 
sorbed by the medium, or their momenta may be shifted, 
it is hard to obtain the reconstruction. Since the meson 
inelastic cross section is supposed to be 2/3 of the baryon 
cross section, the contribution from ITD loss for mesons 
would be more significant than the one for baryons. The 
transverse momentum distributions for p° and to are pre- 
dicted in Fig. [TU 

The 4> production is shown in Fig|TS] The model re- 
sults are much larger than the experimental data. This is 
partially due to the limited hadron list in our treatment. 
Otherwise, the emission at the hypersurface for <j) may 
not be effective. Note the fact that the STAR data [53 
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FIG. 9: The transverse momentum distributions of A and H~ . 
The decay contribution from S~, S° and f2 to A is removed, 
but the one from E° is included. The distribution for H~ is 
inclusive. All the data are from STAR Collaboration [48l.l49l|. 



FIG. 11: The transverse momentum distributions of A, K + 
and Kg. The distributions for K% and "A all" are inclusive, 
and the other decay treatment is the same as in Fie [9] All 
the data are from STAR Collaboration [iEI. lij. |49|. l5o| . 
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FIG. 10: The transverse momentum distributions of A and 
H + . The treatment for decay contribution is the same as in 
FigH The data are from STAR Collaboration |H, Bl- 
are about twice the PHENIX data [53|, further precise 
measurement for <p is required. 

The E*(1385) distribution is shown in Fig. OH The 
model results are much smaller than the experimental 
data, especially in very low pt region. The situation 
here is similar to the case for A. The A-7T or NK pro- 
cess which is excluded in the model may regenerate some 
E*(1385) [H, [H]. As £*(1385) can not decay to NK 
due to its low mass, it may become a one-way valve to 
deliver strange quarks from K to A(S). Another possible 
reason for the E* suppression is that many heavy reso- 
nances which would decay to S*(1385) are not included 
in our calculation. 

The transverse momentum distributions of S*°(1530) 
and f2 are shown in FigQ~7] For fl, the calculated distri- 
bution is close to the data in low px region but deviate 
from the data in high px region. The situation is similar 



FIG. 12: The ratio A/K s as a function of transverse mo- 
mentum. The abundances are inclusive, and no feed-down 
correction is applied. The data are from STAR Collaboration 




to <j). This shows that either more resonances should be 
considered or the scenario of the hypersurface emission 
for these hadrons need to be corrected. From the data, 
it looks that there are two slope parameters for f2, like A 
and S*(1385). As for S*°(1530), the model calculation is 
bad at low px- It might be due to the fact that 5(1950) 
is not counted in our model. S(1950) is a mysterious 
baryon. Its spin, parity and branch ratios are not clear, 
and it is even suggested that there might be more than 
one 5 at that mass |36l |. 

The ratio of £1 to cf> is presented in Fig[TH]as a function 
of transverse momentum. As multi-strange hadrons, the 
jet contribution to <fi and f2 is expected to be small, as we 
have seen in their spectra in Figs [T5l and [T71 Our results 
are similar to the ones calculated in recombination model 
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FIG. 13: The inclusive m T distribution of (if* (892) + 
K*° (892))/2. The dotted line is the model calculation with- 
out considering the in-transport decay (ITD) loss correction. 
The data are from STAR Collaboration JjOj]. 



FIG. 15: The inclusive m T distribution of 0(1020). The x 
error bar is used here to show the rriT bin size. The data are 
from STAR H2] and PHENIX [H Collaborations. 
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FIG. 14: The model predicted inclusive distributions of p' 
and u. 




Some of the momentum integrated ratios are listed in 
Table Mil and compared with experimental data [45], • 
The data from [57| correspond to baryons in the pt region 
of 2.4 — 3.6 GeV and mesons in the pr region of 1.6 — 2.4 
GeV. The influence of hadron decay is relatively small 
in these regions. The strange baryon ratios at colliding 
energy y/s = 200 GeV I57| are very similar to the data 
at energy ^/s = 130 GeV |58l. [H, [6(j. Some meson and 
baryon ratios as functions of transverse momentum are 
shown in Figs fl9l and l20l 

£ and A particles are sensitive to the hadronization 
mechanisms. In some string fragmentation models the 
ratio E°/A is taken to be ~ 0.35, but in thermal sta- 
tistical models it is predicted to be 0.65 — 0.75 due to 
the mass difference. If the hadronization is dominated 
by gluon junction or coalescence, the ratio is about 1. 
In our case, the hadronization mechanism is similar to 



FIG. 16: The inclusive m T distribution of E* ± (1385) + 
E* T (1385). The data are from STAR Collaboration [H|. 

the coalescence or recombination, but the ratio is a little 
larger than 1 [26|, [35| ■ The reason is that S° contains an 
important component Bi (A uc i, s), described in Appendix 
lAl but this term dos not exist for A with isospin zero. 

With the selected parameters, the direct ratio S°/A is 
1.181. Considering the decay contributions from heavy 
baryons or resonances, it decreases to 0.636. The trans- 
verse momentum distributions of direct S° and A, inclu- 
sive E°, and exclusive A without E° decay contribution 
are shown in Fig l21l The transverse momentum depen- 
dence of the ratio E°/A with and without considering the 
decay contribution are shown in Fig. [22j 

Since it is difficult to extract E° through electromag- 
netic decay, E + may be easier to be determined if ir° is 
well measured. Its production is supposed to be similar 
to E°. The expected E + production is shown in Fig. [521 

A(1405) 0(|") and A(1520) 0(§") are regarded as sin- 
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FIG. 17: The inclusive transverse momentum distributions of 
H°(1530) and fi _ + Q + . The data are from STAR Collabora- 
tion mmm. 
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tion \M- 



glet states. It is not easy to determine their wave func- 
tions in quark-diquark model, and it is hard to distin- 
guish them from the states (L = 1,5 = §) and [L = 
1,5 = i). A rough estimation of A(1405) and A(1520) 
production is shown in Fig. [24j where A(1405) and 
A(1520) are, respectively, supposed to be (L = 1,5 = |) 
octet and (L = 1,5 = §) octet. The details on their 
wave function construction is given in Appendix [XJ 

At the end of the discussion on hadron spectra, we 
consider the ratio Cl + /f2 _ . When diquarks are considered 
as constituents of the strongly coupled QGP, the strange 
chemical potential of the system will increase due to the 
negative strange numbers of diquarks us and ds. In our 
calculations, the strange chemical potential /is = 9.35 
MeV is extracted from the strange charge conservation. 
In this case, the strange quark and diauark ss will carry 



Particles 


M odel 


Exp 




7T~/7T + (in) 


1.006 


0.984 ± 0.004 ± 0.057 


[45] 


7T~/7T + (ex) 


0.999 


1.01 ±0.02 [57] 




K-/K+ 


0.915 


0.933 ± 0.007 ± 0.054 


[45] 






0.96 ±0.03 [57] 


PAP(ex) 


0.748 


A 701 _L A fl1 1 _L A Ai^O 

0.731 ± 0.011 ± 0.062 


145] 


71 / Tl I 111 1 


0.783 


n 747 _|_ n flfl7 + o 046 


1451 






a 77 _i_ a ac rc7i 


X / A 


U.o44 


n 70 -i- n no/i r^7l 

U. i Z It 






0.907 


0.82 ± 0.05 [57] 




fr/rr 


1.011 


1.01 ±0.08 [57] 




K+/TT + 


0.203 


0.171 ±0.001 ±0.010 


[45] 


t\ 1 71 




n i fio 4- n nm -4- n ni n 






0.0645 


0.064 ± 0.001 ± 0.003 


[45] 


p/7r+(in) 


0.102 


0.099 ± 0.001 ± 0.006 


[45] 




0.0479 


0.047 ± 0.001 ± 0.002 


[45] 


p/7T"(in) 


0.0797 


0.075 ± 0.001 ± 0.004 


[45] 



TABLE III: Some momentum integrated hadron ratios com- 
pared with RHIC data [45l [57j . The statistical errors in the 
model calculations are up to 5%. 
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FIG. 19: The inclusive n /7T and exclusive p/p ratios 
as functions of transverse momentum. The data are from 
PHENIX Collaboration El. 



negative chemical potentials, fi s — —0.35 MeV and/i ss = 
—0.7 MeV [35[. The string fragmentation model [6lj with 
diquarks gives similar result. However, the ratio in the 

~" < 1 sin. 

7T+ 



calculation without diquarks is Q /fi 

As shown in Table [m] the ratio TV /Q,~ = 1.011 
is consistent with the STAR measurement [57| . While 
the diquarks are supposed in thermal equilibrium in our 
model, and the current data are still with large statis- 
tical and systematic errors, the qualitative conclusion of 
f2 + /Sl~ > 1 is supposed not to be changed. 

As pointed above, the thermalization of final state 
hadrons are omitted in all of the calculations. Consider- 
ing the interactions which are not inserted in our treat- 
ments, some of the dominant hadrons would be partially 
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FIG. 22: The ratio E°/A as a function of transverse momen- 
tum, with and without considering decay contribution. 
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FIG. 21: The model predicted transverse momentum distri- 
butions of direct E° and A, inclusive E°, and exclusive A. The 
decay contribution from E°, E~, H° and fi to the exclusive A 
are removed. 



thermalized effectively. Comparing with the RHIC data, 
7T, K, p and A may be well thermalized, while the heavy 
baryons may be not. 



V. SUMMARY 

We have developed a phenomenological model to de- 
scribe hadron production in relativistic heavy ion col- 
lisions. In our model, the QGP evolution and hadron 
transport are coupled to each other in a self-consistent 
way. The strongly coupled QGP, including not only glu- 
ons and quarks but also their resonant states like di- 
quarks, is described by ideal hydrodynamics, and the 
hadrons emitted from the QGP satisfy transport equa- 
tions with loss and gain terms. The hadron decay is ef- 



FIG. 23: The expected transverse momentum distributions of 
direct E + and E~. No feed-down correction is applied. 



fectivcly reflected in the medium absorption and hadron 
spectra. Due to the feedback of the hadron emission, the 
hadronization hypersurface of the parton fluid shrinks 
with a constant velocity. The splitting of the fluid due to 
its fast expansion in the later stage and the hadron emis- 
sion is also effectively considered, by introducing a ratio 
of the volume of the QGP droplets to the whole volume 
of the system. We applied our model to calculate the 
final state hadron spectra in heavy ion collisions. From 
the comparison with the RHIC data, most of the hadron 
spectra and hadron ratios can be reasonably described in 
our model. 

The model deserves further improvement in future 
studies. Unlike quarks and gluons, diquarks in QGP are 
resonant states and they should behave differently. As a 
consequence, the assumption that the baryon cross sec- 
tion is 1.5 times the meson cross section, which arises 
from the constituent quark model, should be reconsid- 
ered carefully. Since mesons are described by transport 
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FIG. 24: The model calculated transverse momentum dis- 
tributions of A(1405) and A(1520), with the wave functions 
listed in Appendix [A] No feed-down correction is applied. 



equations in our model, the assumption of diquark ther- 
malization is in principle not self-consistent. Different 
types of hadrons will have different thermalization time. 
Light hadrons are easy but heavy hadrons are difficult to 
get thermalized. In our current treatment, the relaxation 
time in hadron transport equations is taken to be infin- 
ity, and therefore, no hadron can be thermalized in the 
whole evolution. The treatment of the QGP fluid split- 
ting is also too simple in our model. When the model is 
improved, it can be used to calculate more hadron dis- 
tributions like elliptic flow 63, 64], which is considered to 
be sensitive to the thermalization of the parton system 
in the early stage of heavy ion collisions. 

Acknowledgements: We are grateful to Wojciech Bro- 
niowski, Huanzhong Huang, Rudolph C. Hwa, Yuxin Liu, 
Sevil Salur, Qubing Xie, Nu Xu and Weining Zhang for 
useful discussions. The work is supported in part by the 
grants No. NSFC10547001 and 10425810. 



APPENDIX A: REWRITING THE BARYON 
WAVE FUNCTIONS 

The SU(6) baryon wave functions can be rewritten in 
the quark-diquark model [H, [66| . The (^) + and (|) + 
s-wave intermediate baryon states are introduced as 



B h+i {A,q) = 
D 



1 

71 

Sq\ 

AW 
1 



V2A+q l - A q ] 



+ ±{A 7 q) = A, ,r +V2A q 1i 



(Al) 



where S and A represent a scalar and an axial-vector 
diquark state, respectively. The final baryon states can 



be expressed in terms of certain combinations of these 
intermediate baryon states [26l [35| . The wave functions 
are expressed as 
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4 Bx(A us ,d) - d -Bx(A ds ,u) 



+J\B i (S us ,d)-J±Bi(S ds ,u) 



+Bi(S ud ,s) 
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2 Bi(A ds ,s) 
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(A2) 
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for the SU(3) octet baryons, and 

Bs (A uu , d) 
B3(A dd ,u) 



|A++) 
|A+) 

|A°) 
|A") 
|E*+) 

|E*°) 

|E*") 

|S*°) 

| s*-) 

10") 



-Ba {A uu , . 
1 r 



V3 
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1 r 

71 
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71 
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^(A^s) 



V2B3(A ud ,u) 
V2Bs_(A ud ,d) 



B3(A uu ,s) + V2B3(A us ,u) 
B 3 .(A us ,d) + B^(A ds ,u) + B i (A ud ,s) 
B i (A dd ,s) + V2B i (A ds ,d) 
B3_(A ss ,u) + V2B3(A us ,s) 
B3{A ss .d) + V2B3{A ds ,s) 



(A3) 



for the decuplet baryons. From these expressions, we can 
get the factor Cj in Eqs. (|2"T)) and (|3"Tj) to complete the 
production rate. 

There might be other orthogonal baryon states. From 
the Tables of Particle Physics [H, 7V(1440) ±(± + ) and 

N(\7\0) \{\ + ) are conjectured to be the orthogonal 
states of p and n. Other possible orthogonal states 
are A(1600) 0(± + ), A(1810) 0(± + ), E(1660) 1(± + ) and 
E(1880) l(i + ) of A and E, AT(1720) |(| + ) f A+ and A , 
A(1890) 0(§ + ) and E(1840) 1(| + ) of E* (1385), and some 



unknown orthogonal states of 5 and 3* (1530). There are 
24 such particles which belong to two {\) + octets and a 
(§)+ octet. 

The wave functions of the orthogonal baryon states in 
(|) + decuplet can be determined. For example, 



|A 2 )^ } 3 = -j= [B3_{A us ,d) - Bs 2 (A ds ,u) 



(A4) 



|S^) { „ 8 _ } 3 = \B3(A us ,d) + B3(A ds ,u) - 2B3(A ud ,s) 
2 v6 L 2 2 2 

There are two orthogonal baryon states in (|) + octet. 
Their wave functions are not unique. Considering that 
the first octet state is the mixed state of Aq and Sq, a 
most possible assumption is that the second octet state 
should be the corresponding mixed state, and the third 
one contains no mixing. Therefore, we have 



1 
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V3 



Bi(A uu ,d) - 
^ Bi(S ud ,u) 



Bi_{A ud ,u) (A5) 



Bx(A uu ,d)+V2B l (A ud ,'< 



where, N% corresponds to A^(1440) + and to 
7V(1710) + . The other orthogonal baryon wave functions 
with positive parity are similar. 

The wave functions of some of the and (|)~ res- 
onances are not clear. It is hard to distinguish a {\)~ 
or a (|) _ state from the two candidates (L = 1, S = 5) 
and (L — 1, S = §). For simplicity all (i)~ baryons, ex- 
cept A, are regarded as S = \ octet states and all (|)~ 
baryons are regarded as S — § octet or decuplet states. 

A(1405) 0(i~) and A(1520) 0(§~) are regarded as sin- 
glet states [3g | . The quark-diquark wave functions of 
them are unknown. Approximately, Their wave functions 
are written as S = \ octet and S = | octet. 



APPENDIX B: EFFECTIVE LAGRANGIAN AND 
MATRIX ELEMENTS OF PRODUCTION CROSS 
SECTION 

In calculating the hadron production rate, only low- 
est order interactions are considered. As there are not 
enough data to determine the coupling modes and con- 
stants, we choose only some of the most simple coupling 
channels. To further reduce the number of parameters, 
we use two universal coupling constants gM for meson 
interactions with two quarks and g B for baryon interac- 
tions with a quark and a diquark. The meson sector of 
the effective Lagrangian density is listed in Tab II VI for 
different meson. The influence of P-wave and D-wave is 
neglected. 



Meson 



\M S 



^ [m 2 M - (m q i - m q 2) 2 ] 
a -f- [m 2 M - (m ql +m q 2) 2 ] 

a\l [ m M-( m gl ~m q2 ) 2 



o~ ig M qijsq24> 

2 + ig M <Zi7m7^2K m1 , %^ \m 2 M - (m q i + m q2 ) 2 ] 



9M [ m2 M -( m ql+ m q2) 2 ] 



I m M 
I m M 



+ 2j 

+ 2] 



TABLE IV: The effective interactions and the corresponding 
scattering matrix elements |M sp in| 2 = |M| 2 /3 for different 
mesons. 

For the baryon sector of the effective Lagrangian den- 
sity, we take 



L q sB = ifjB BqS, 
L q ABi = i9B Bj^qA^, 
\M\ 2 S = g 2 B [(m B + m q f - m%] , 



(Bl) 



M\\ = 9 -f[{ mB+mq f 



(m B - m q ) 2 



m 



for the baryons with J 1 



i + , L = and S 



Note 



that the scattering matrix elements here are equivalent 
to the expressions in [35j. 
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For J p 

LqSB — 



|67| . L = 1 and S = |, we use 



qABl 



\M\i = 



m B +m q 
igB 



(B2) 



|M| 



2m a 

9b [( m B - m q ) 2 - nig] , 

9 b [i \2 2 

— [(m B - m g J - m A 



(rn B + m q ) 2 



with the definition — —i[j f j - ,j u ]/2. The expression 
of L q gs here gives the same result as the direct coupling 
(6^ | ■ For J p = | , L = 1 and S = |, the scattering 
matrix elements are supposed to be the same. 
We choose 



L q AB = igB B^qAft, 



M\ 2 A = ff [{m B - m q f - m\ ] 
lo 



(rn 



,2\2 



"B A 



(B3) 



for J 1 



! + 



, [69|], L = and 5 = |, and 

«5s 



mg 

9l 



B^qd^S, 



\M\% = f[(m B -m q f-m s ] 



m 2 f 



- 4 



B" L S 



(B4) 



couplings, the matrix element is replaced by 



\M\\ = f| [(m B +m q f 



(m| 



m 2 q f 



l B" l A 



(B5) 



and the expression for J p — | , L = 1 and S = | is the 
same. 

We employ 



L q ab — — — Ba V qduA Vl 
m A 



M\\ 



|| [(m B + m,) 2 - ™ 2 



2\2 



for J 1 



"B'^A 
^,2 



5.6 



1 and S 1 



(B6) 



|, where the expression of 



the spin-| projector can be found in [70|, and 



JqSB 



— Y Buvqdud u S, 

1% 



mi 



= ^[(m B +m q ) 2 ~m 2 ] 



(m 



1% - m 2 q ) 2 



B" b S 



(B7) 



for J p = 



L = 1 and 5 = \. For the axial- vector 



for J p = | + ,i = 2and5=i. 
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